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We describe and systematically analyze new implementations of the Projection 
Operator Imaginary Time Spectral Evolution (POITSE) method for the Monte 
Carlo evaluation of excited state energies. The POITSE method involves the com- 
putation of a correlation function in imaginary time. Decay of this function con- 
tains information about excitation energies, which can be extracted by a spectral 
transform. By incorporating branching processes in the Monte Carlo propaga- 
tion, we compute these correlation functions with significantly reduced statistical 
noise. Our approach allows for the stable evaluation of small energy differences in 
situations where the previous POITSE implementation was limited by this noise. 



1 Introduction 

The Projection Operator Imaginary Time Spectral Evolution (POITSE) 
method has allowed calculation of excited states_to be made with diffusion 
Monte Carlo (DMC) without nodal constraints III The main requirement is 
that a reasonable ground state wave function be available, which can be ob- 
tained from well-established ground state methods such as DMC. The excited 
states are then accessed via projector operators, whose evolution in imaginary 
time contain information on excited state energies. In the POITSE method a 
correlation function of the projection operators is evaluated by Monte Carlo 
techniques, and then subsequently inverted to obtain spectral functions whose 
peak positions correspond to excited state energies. This inversion requires 
an inverse Laplace transform, a notoriously ill-conditioned numerical proce- 
dure. In the applications of POITSE made to datejalffl this inversion has been 
performed with the Maximum Entropy Method (MEM).! POITSE has con- 
siderable power in allowing analysis of excited states without imposing nodal 
restrictions. It is particularly useful when some physical insight about the 
nature of the desired excited states is available. This information can be used 
to tailor suitable projectors to obtain maximum overlap with the eigenstates 
of interest. This has been demonstrated recently with permutation symmc- 
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try tunneling excitations.u In general, the viability and power of the method 
has now been shown for a range of model systems involving atomic motionsB 
It has been applied to several physical examples of cluster excitations which 
cannot be^ addressed by basis set methods, including up to 15-dimensional 
problems!! To our knowledge, the method has not yet been systematically 
applied to fermion problems, although there is no intrinsic impediment to 
this. 

In this paper we analyze the efficiency and accuracy of the POITSE al- 
gorithm for various different implementations of the DMC component of the 
method. We present a modification of the algorithm that allows the calcu- 
lation of small energy differences with reduced statistical noise. In Sec. ||, 
we briefly review the POITSE general formalism and explain in detail the 
different numerical implementations. Sec. || illustrates the different imple- 
mentations with two applications: the one-dimensional problem of the am- 
monia inversion mode and the six-dimensional van der Waals vibration of the 
4 He-benzene dimer. 

2 Computational Methodology 

The general POITSE method involves the Monte Carlo evaluation of an 
imaginary time (r = it) correlation function k(t), and then a subsequent 
inverse Laplace transform of this correlation function using the Maximum 
Entropy Method. With an appropriately chosen correlation function, the in- 
verse Laplace transform provides a spectral function whose peak positions 
correspond to ^excitation energies. The basic theorjJd and its application to 
model systemsB'ETQ have previously been described in detail, and thus we will 
only present a brief summary of the relevant formalism. 

2. 1 Theory 

The primary quantity of interest in POITSE is the spectral function k(E), 
k(E) = \(MA\<f>n)\ 2 S(E -E n + E Q ), (1) 

n 

where {| <f> n }} and {E n } are a complete set of energy eigenkets and eigenener- 
gies for the Hamiltonian H , and A is an operator chosen to connect |(?!>o) at 
least approximately to the particular excited state(s) of interest \4> n ). Taking 
the Laplace transform of Eq. (|]J) , one can obtain the imaginary time correla- 
tion function k(r), in atomic units (h = 1), as 

k(r)^(MAe^ A - E ^A^\^) (2) 
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The POITSE approach consists of evaluating k(t) by a Monte Carlo algo- 
rithm, then taking its inverse Laplace transform to obtain the spectral func- 
tion k(E). 

In most situations, however, the ground state \(f>o) is not known exactly. 
In practice, one typically employs a trial function I^t) and reference energy 
E re f which approximate as closely as possible \4>q) and Eq, respectively. Use 
of a reference energy not equal to the exact ground state energy modifies the 
decay rate of all terms in Eq. (||) by a constant factor of E re f — Eq. This 
results in a systematic bias in the excitation energies of Eq. ([!]) , independent 
of the usual finite time step bias due to the DMC evaluation of Eq. (||) . This 
bias from E re f is also independent of whether the true ground state |^>o) is 
used. n 

It has been shown earliertl that such systematic bias can be eliminated by 
introducing the normalization factor 

(* T |e-^- B "^ r |* T ). (4) 

The removal of bias due to E re t can be seen from the following arguments. 
First, replacing \4> ),Eq in Eq. (||) with |\I/t), -Eve/, respectively, and dividing 
by the additional normalization factor of Eq. , leads to the modified decay 
function 

\^t) is then expanded in eigenstates of H to yield0@ 

E„|(*T|i|0„)| 2 e- (iJ "-^ )r 



k(t) 



(G) 



c 2 e -(E m -E„ f ) T 

where c m = {^T\4>m)- The numerator and denominator of Eq. (|^) may then 
be multiplied by e ( - E °~ Er ^ T /c^ to obtain 
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oc £ |(*T|i|^„)| 2 e- (E »- So)T + 0{x). 

n 

Here, the prefactor of Eq. ([?]) was expanded in a power series in x, where 
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When I^t) = \(f>o), we see that Eq. (||) is identically equal to Eq. (||), and the 
effects of using a reference energy other than the true ground state energy are 
completely eliminated. Additive errors of O(x) and higher are present when 
an approximate ground state is used. Note that since the series expansion of 

Eq. (0) is only convergent for c > ^J\, this does require that a reasonable 

approximation to the ground state be available. The higher order terms 0(x) 
contribute to the spectral function k(E) in an additive manner. Consequently, 
they do not affect the positions of the relevant spectral features of interest, i.e. 
the dominant leading terms of Eq. (Js[) . In practice, for a xeasonable choice of 
\^t) these additional terms have highly reduced weight.Eltl To leading order 
therefore, the renormalized decay Eq. (||) exhibits the time dependence of 
Eq. (|J), independent of the reference energy E re f. Consequently E re j may 
be arbitrarily chosen and varied. The usefulness of this will become more 
apparent below, in our discussion of numerical implementation. 

The numerical inversion of H(t) to obtain k(E) is an ill-conditioned prob- 
lem, especially when Monte Carlo noise is non-negligible and/or when the 
spectral function k(E) contains multiple overlapping peaks of comparable in- 
tensity. Thus a judicious choice of the operator A is necessary to ensure that 
the time-dependence of k(t) is dominated by only one or a few well-separated 
energy differences. The imjerse Laplace transform of k(t) is performed using 
the Bryan implementations of the maximum entropy method. Our use of 
this approach for the. imgrsion of k(t) is identical to that employed in pre- 
vious POITSE work.tffltm We will discuss choices for A specific to particular 
systems of study in Sec. ||. 



2.2 Numerical Implementation 

The correlation function q£ Eq. (^|) may be rewritten in a form amenable to 
Monte Carlo evaluation aslil 

K(T)=£ ,A.(Rf)i(RrM<», i m 

(t) 

where is a guided random walk j in multidimensional configuration space, 
discretized in time steps of size At (a DMC "walker"), and 

W (R< r) ) = nexp{-[£L(R^ mAT) ) - E ref ]Ar}, (11) 

m 

ElCR^) = * T 1 (R( T) )i?vl/ T (RM). (12) 
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The quantities w(Rj ) and _EY(R^ ) are the usual DMC cumulative weight 
and local energy, respectively.0 The evaluation of Eq. ( |l0| ) begins with a varia- 
tional Monte Carlo (VMC) walk in which an initial starting ensemble of walk- 
ers distributed according to ^j.(R) is generated using a simple Metropolis 
method.El The starting VMC ensemble is subsequently propagated in imagi- 
nary time by a DMC sidewalk, during which Eq. ( |l0| ) is sampled. Since the 
maximum entropy analysis requires independent samples of H(t), the start- 
ing configuration for each DMC sidewalk is taken from the VMC walk every 
100 — 200 VMC steps apart, to minimize correlations between successive side- 
walks. The set of re(r)'s evaluated in this manner serve as input for the 
inverse Laplace transform via MEM. Typically 100 — 500 independent decays 
are required to produce a converged spectrum k(E). 

In the original implementation of Blume et aZ.,0 the DMC weights w(R^) 
take on a continuous range of values, and walkers are not destroyed or du- 
plicated. We refer to this approach here as DMC with pure weights. This is 
the preferable implementation in an ideal situation where high-quality trial 
functions are available. However, for reasonably complex systems this is often 
not the case. In addition, it has been shown that DMC with pure weights is 
unstable for long propagation timesB Therefore, as we demonstrate in Sec. ||, 
a DMC sidewalk that uses pure weights may sometimes be impractical in 
situations involving small energy differences. 

A common solution to the problems associated with pure weights is to 
introduce branching. The simplest branching scheme rounds the walker weight 
at every step of the walk to an integer rij = int[iu(Rj- ) + £], where £ is an 
uniformly distributed random number on [0, 1). A walker R^ is destroyed 
for rij — 0; otherwise, rij copies of walker R^ 1 "'' are propagated independently 

in the next DMC move. In this case, the weights ?i)(Rj ) take on only integer 
values, and Eq. ( |Io|) becomes 

~k(t) = — £it(Bf )i(RM), (13) 

3' 

where the index j denotes the parent walker at initial time r = from which 
walker j' at time r descended, and the instantaneous ensemble size n w fluctu- 
ates with time. We refer to this approach here as DMC with pure branching. 
While the pure branching method is formally correct on average and is much 
more stable numerically, the integer_rounding of walker weights can neverthe- 
less lead to greater statistical noiseJj 

To minimize this noise, one can employ a hybrid approach where each 
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weight w(R^ T ') is allowed to vary continuously, and a walker is only destroyed 
or duplicated when its weight exceeds some predetermined bounds. In such 
a situation, it is important that the branching procedure does not artificially 
alter the ensemble sum of weights W to t = J2j w(R^). A combined weighting 
and branching scheme will in general exhibit less statistical noise than a pure 
branching scheme. In some cases the noise reduction can be significant. Our 
implementation of branching is similar to that outlined in Ref. About 
every 20 — 50 DMC time steps, the ensemble is checked for walkers whose 
weight exceeds the empirically determined bounds w m i n and w max . A walker 
R^ T) with weight w(Rj) > w rnax is split into nj — int[w(R^ r ' ) ) +£] walkers, 
each with weight w(R.j)/nj. A walker Ry with weight u>(Rj r ' ) ) < w m i n is 

either a) killed with probability 1 — to(R^), otherwise b) kept with its weight 
set to unity. The bounds w m i„ and w max are chosen to give a stable DMC 
walk with respect to the ensemble size and W to t- 

As discussed previously, incorporation of the normalization factor of 
Eq. (Q) into k(t) results in a decay independent of the reference energy E re f. 
Therefore we are free to choose and vary E re f based on considerations of nu- 
merical stability. A common choice of E re f is the variational energy of the 
trial function, E re t = (^>t\H \^t) / (^t\^t), which may be obtained from 
a separate VMC calculation. One may also choose the ground state energy 
E^ = Eq, which is readily obtained from standard ground state DMC meth- 
ods. In our implementation, we begin with an initial choice of E re f and update 
E^ continuously during the course of the DMC walk according to 



Kef ~ ^ref + m 



vW r+Ar) ) 



(14) 



where r\ is an empirical update parameter chosen to be as small as possible 
to avoid biasing the results, typically ij/At = 0.01 — 0.3. The effect of this 
updating procedure for E re f is to keep the average walker weight close to 
unity, thus preventing the ensemble size and sum of weights from diverging 
off to infinity or zero. The combination of these various mechanisms serve 
to ensure a stable DMC walk for long times, thus allowing the evaluation of 
small energy differences E n — Eq. In the examples presented in Sec. ^ we will 
compare the effects of the various DMC schemes described here. 

A final note in the implementation concerns the statistical errors in the 
excited state energy differences E n — Eq- The MEM inversion of k(t) gives 
the spectral function k(E), whose peak positions correspond to excited state 
energy differences. There is no general approach to assign error bars in the 
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mean peak position^ and thus we only report energies to the last significant 
figure. We determine empirically the position of this last significant figure by 
examining the convergence of E n — Eq with respect to the number of decays 
k(t) used as input for the MEM inversion. Because multiple projectors are 
usually sampled from the same DMC sidewalk, the relative differences between 
excited states are expected to be very accurate. 



3 Examples 

3.1 NH3 inversion 

The first application we discuss here is a POITSE study of the ammonia 
inversion mode. Freezing all other internal degrees of freedom, the Schrodinger 
equation for this mode alone is a one-dimensional problem which can be solved 
exactly by a straightforward Discrete Variable Representation-Finite Basis 
Representation (DVR-FBR) calculation. EJ The Hamiltonian is given by 

k =- n i^ +v{hl (15) 

where h is the distance between the nitrogen atom and the hydrogen plane, 
fi is the effective mass for the mode and V(h) is the double-well inversion 
potential for tunneling across the hydrogen plane. We use one of the poten- 
tial forms ("Case b") proposed by Nino et al£3 which leads to a tunneling 
splitting of 1.43 cm -1 for the lowest tunneling pair, and 64.5 cm^ 1 for the 
next lowest tunneling pair. The corresponding DVR-FBR energy levels are 
listed in Table [l] as benchmarks for the POITSE results. 

A double well study was previously made in Ref. ^ to demonstrate the 
effectiveness of the POITSE method for model systems. However, in that 
example, the energy differences involved were much larger than those arising 
in the NH3 inversion problem which we discuss here. While the inversion fre- 
quency is high (993 cm" 1 ), the POITSE method allows the computation of an 
energy difference which is three orders of magnitude smaller. We compare here 
two different DMC implementations, namely pure weights and pure branch- 
ing, and demonstrate the limitations associated with the former approach for 



Table 1 . Lowest four energy levels (in cm ) for the inversion mode of NH3 relative to the 
ground state energy, which is 553.11 cm -1 above the potential minimum. 





Eq Ei E 2 E 3 


DVR-FBR 


0.00 1.43 961.40 1025.93 
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computing small energy differences. In order to make such a comparison of 
implementations, it is convenient and indeed preferable to use a system for 
which exact wave functions can be found. 

The trial function vf^/i) used in the Monte Carlo evaluation of Eq. (|J) 
was initially fit to the DVR-FBR ground state cigcnfunction $o(h), and then 
further optimizfid.Jai.VMC. While numerous sophisticated VMC optimization 
schemes exist Jl3'EJ'[L3 for a simple one-dimensional problem we found it suffi- 
cient to manually vary the trial function parameters to minimize the ground 
state energy and its variance. We use the analytical form 

* T (h) = exp[a e bo(h - co)2 + a e bo(h+Co)2 + d e e ° hi }, (16) 

where do, bo, cq, do, and eo are parameters listed in Table ^. The corresponding 
VMC energy is 561.4(3) cm -1 , which is less than 2% above the exact ground 
state value obtained from DVR-FBR. 

Since the first excited state of a double well potential is the lowest 

antisymmetric state, the projector A — h was previously usedo to access this 
level. In obtaining higher- lying states, choosing A to be an integer power of h 
led to a k(t) consisting of a superposition of multiple exponential decays. For 
instance, a choice of A^T(h) = h 2 ^x{h) resulted in non-negligible overlap 
with multiple excited levels. Thus an accurate Laplace inversion of the corre- 
sponding k(t) was more difficult, due to the multiple decay contributions of 
these states. We use here instead more effective projectors given by the ratio 
of the eigenfunctions 

where $ n (/i) is an excited state eigenfunction obtained from a DVR-FBR 
calculation. Clearly if the eigenfunctions are numerically exact, this results in 
an exact projector. Such projectors have also been shown to be useful n when 
only symmetry properties of the eigenfunctions are well characterized. □ The 
following analytical expressions were fitted to the DVR-FBR eigenfunctions 
for the lowest three excited states (n — 1 — 3): 

$ x (ft) = e bl(h - Cl)2 - eM h+c ^ 2 (18) 
* 2 (/0 = a 2 [(h - f 2 )e b ^ h - c ^ 2 -{h + f 2 )eM h+c ^ 2 } + d 2 e e ^ (19) 
$ 3 (/0 - a 3 [(h - h)e H{h - C3)2 + (h + h)e b ^ h+c ^] + d z he^ h \ (20) 

The fit parameters are given in Table |[ We emphasize that we are using 
this example of ammonia inversion to demonstrate and compare the relative 
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Table 2. Fit parameters (in atomic units) for NH3 eigenfunctions obtained from DVR— FBR. 



7) 


a n 


K 




d n 




fn 





17.0 


-1.095 


0.829 


1. 


-0.054 




1 




-10.886 


0.681 








2 


0.785 


-11.072 


0.674 


-0.0664 


-0.082 


0.580 


3 


0.768 


-12.025 


0.720 


-0.314 


-1.325 


0.447 



efficiency of two different and alternative implementations of the POITSE al- 
gorithm. In particular, we shall compare the extent of noise and time step 
bias of the two different approaches to the DMC evaluation of Eq. (J5J) . Our 
aim here is not the establish the generality of the method, or its accuracy far 
a double well problem, both of which have been addressed in earlier workJilu 
Instead, we are interested in assessing the relative efficiency of different im- 
plementations, and thus it is preferable here to use projectors which are as 
exact as possible. 

Since the lowest tunneling splitting E\ — Eq is small, the correspond- 
ing decay k(t) is slow and requires a long DMC propagation. Fig. |l|a shows 
four typical k(t)'s computed using the original POITSE implementation in- 
volving DMC with pure weights. These decays become extremely noisy as 
the time r increases. The ensemble Jpcal energy (El) also exhibits such be- 
havior. This problem is well-knownLl and arises from the fact that for long 




"0 0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 



i[Hartree 1 ] x io 5 x[Hartree 1 ] x io 5 

Figure 1. Typical correlation functions k(t) for NH3 using the projector ^i/ffj". The left 
plot (a) corresponds k(t) evaluated using DMC with pure weights, while the decay curves 
in the right plot (b) are obtained using DMC with pure branching. 
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or even moderate DMC propagation times, the Monte Carlo ensemble aver- 
ages are dominated by only a few walkers carrying high relative weights. In 
comparison, Fig. |l|b shows typical k(t)'s obtained from an implementation 
using DMC sidewalks with pure branching, where walkers are replicated or 
destroyed at each time step based on integer rounding of their weights as dis- 



cussed in Sec. 2.2. In both calculations, 2000 walkers were propagated using 
a time step At of 5 Hartree -1 . Clearly there is far less noise at longer times 
in the pure branching implementation, and thus such an approach is more 
suitable for the evaluation of small energy differences. Using the pure branch- 
ing scheme, the Laplace inversion of 600 decays computed up to a final time 
Tf of 250000 Hartree -1 results in a single peak at 1.39 cm -1 , in reasonable 
agreement with the DVR-FBR value. 

The evaluation of the larger energy differences E% — E and E 3 — E are 
manageable using both DMC implementations, because the lengths of the 
corresponding decays are much shorter than for the lowest energy difference 
Ei — Eq. The use of the projector given in Eq. ( |l7| ) facilitates the Laplace 
inversion, since each choice of A n results in a k(t) consisting of only one 
exponential decay. In these calculations, 1000 decays are used as input for 
the MEM inversion, with each decay computed using an ensemble of 2000 
DMC walkers propagated to a final time tj of 1500 Hartree -1 . The number 
of decays required for a converged k(E) depends on the energy difference of 
interest and on the time step At. In general, for larger time steps, DMC with 
pure weights requires more sampling to produce fully converged results. 

Since DMC methods are subject to a systematic time step bias, we per- 
form a comparative study of the two implementations and their time step 
dependence. For the computation of the lowest energy difference E\ — Eq 
using DMC with pure branching, we find a time step of 5 Hartree -1 to be 
sufficiently small to give an accurate result within statistical error. However, 
the time step dependence of higher energy differences is not necessarily the 
same as that for E± — Eq. Fig. || presents the time step dependence for the 
calculation of E 2 — E and E 3 — E , using both DMC with pure weights (solid 
circles) and DMC with pure branching (open diamonds) . It is evident that for 
both DMC implementations, the higher energy differences are more sensitive 
to time step bias than the lowest energy difference, Ei—E . Thus, in order to 
extract the correct energies in the higher energy range, either a smaller time 
step would need to be used, or an extrapolation to At = would need to be 
performed. 

With this simple example, we have shown that two different POITSE im- 
plementations, namely DMC with pure weights and DMC with pure branch- 
ing, lead to the same results. We have also presented a systematic study of the 
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5 7.5 10 

At [Hartree~ 1 ] 



5 7.5 10 

At [Hartree -1 ] 



Figure 2. Time step dependence for the energy differences E2 — Eo (left) and E3 — Eq 
(right) of NH3 inversion mode. The dashed lines correspond to the exact DVR-FBR values. 
Energies obtained from DMC with pure weights are marked with filled circles, and energies 
obtained from DMC with pure branching arc marked with open diamonds. 



convergence behavior for these two different approaches, and compared with 
the exact solution obtained from DVR-FBR calculations. For the evaluation 
of small energy differences, we conclude that a pure branching DMC sidewalk 
is considerably more efficient than using DMC with pure weights. 

3.2 4 He-benzene dimer 

We now demonstrate the use of the POITSE approach for the computation of 
excited vibrational energies of the 4 He-benzene dimer. We treat the benzene 
as a rigid molecule, and for simplicity we also neglect the rotational kinetic 
energy of the benzene, i. e. the rotation of benzene relative to helium. In the 
space-fixed frame, the resulting Hamiltonian is 

H = -^Vl-fvl + V(r), (21) 
zmn 2m 

where tuq is the benzene mass, m is the helium mass, Vq is the Laplacian with 
respect to the benzene center-of-mass position ro, Vj. is the Laplacian with 
respect to the helium position r^, and V(r) is the 4 He-benzene interaction 
potential. The latter depends only on the relative coordinate vector r = 
Yk — ro. Thc_potcntial is an analytical fhx£l to ab initio MP2 calculations of 
Hobza et a/.JHI and possesses two equivalent global minima of —66.01 cm -1 
along the six-fold Cg-axis, situated at 3.27 A above and below the benzene 
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plane. While in principle one could transform the Hamiltonian to the center- 
of-mass frame to yield a three-dimensional problem, as would typically be 
done in a basis set calculation, sampling the transformed kinetic energy terms 
becomes more complicated in DMC as additional particles are added, and thus 
it is technically simpler for us to work with the six-dimensional Hamiltonian 
as written in Eq. 

The trial function (r) is the product of an anisotropic Gaussian binding 
factor centered on the benzene center-of-mass, and an atom-atom repulsive 
factor, 

* T (r) = e -«(* 2 +!/ 2 )-c Z 2 JJ e *-(r-) JJe^M, ( 2 2) 

a 

where we use for the binding parameters (in atomic units) a = 0.05, c = 0.06. 
The product over a and f3 runs over the carbon atoms and hydrogen atoms, 
respectively. The atom-atom terms t a (r a ) and tp(rp) are functions of 4 He- 
carbon and 4 He-hydrogen distances r a and rp respectively, and their analyt- 
ical forms are chosen to cancel out the leading singularities in the atom- 
atom potential energy terms.ta In this study we use t Q (r Q ) = — c a r~ 6 , 
tpirp) = —cpr^ 5 , with the parameters (in atomic units) c a — 6000, eg = 8000. 
The trial function of Eq. (|2^) possesses the same D^h symmetry as the 4 He- 
benzene potential. 

A ground state DMC calculation using the trial function and potential 
discussed above gives a ground state energy Eq = —21.61(2) cm -1 , which 
corresponds to about 32% of the global energy minimum of the 4 He-benzene 
potential-. Such a high zero-point energy is typical of helium van der Waals 
systemsJiS and underscores the need for a fully quantum mechanical treatment 
of the van der Waals degrees of freedom. 

We choose the excitation operators A^ r ' based on symmetry considera- 
tions, where the superscript T denotes an irreducible representation of the 
Dsh point group. Since the trial function ^t(y) transforms as the totally 
symmetric representation A\ g , for a given A^ r \ the integral (^tI^ 1 ^ \<f> n ) m 
Eq. (^|) is only nonzero for states \(f>n) which transform as T. Thus an ap- 
propriate choice of A^ r ' will, by symmetry, significantly reduce the number 
of terms in the summation of Eq. (^), leaving only decay terms whose char- 
acteristic decay times are presumably more well-separated, and thus easier 
to resolve. The various choices of the operators A^ we use here are listed 
in Table ^, where A^T) is defined with respect to the benzene principal axis 
frame centered on the benzene center-of-mass. In this coordinate system, the 
x-axis is perpendicular to the benzene C-C bond, the y-axis lies along the 
benzene C-H bond, and the z-axis is perpendicular to the benzene plane. 
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Table 3. Operators A^ T ' and the resulting energies E — Eq (in cm" 1 ) for 4 He-benzene van 
der Waals excitations. For the two-dimensional irreducible representations, the two pro- 
jectors listed give degenerate energies. The three rightmost columns list energies obtained 
from hybrid branching/weighting (B/W), pure weighting (PW), and pure branching (PB). 



r 


A?) 


E-E Q 


B/W 


PW 


PB 


Elg 


xz, yz 


6.39 


6.39 


6.39 


Eiu 


x,y 


7.04 


6.97 


7.04 


A 2u 


z(x 2 + y 2 ) 


7.76 


7.64 


7.86 


Alg 


x 2 + y 2 


8.44 


8.54 


8.44 


Eiu 


z(x 2 — y 2 ),xyz 


9.41 


9.36 


9.48 


Elg 


x 2 - y 2 ,xy 


9.96 


9.84 


10.01 


B 2 u 


x 3 — 3a;?; 2 


11.22 


11.34 


11.19 


Big 


z(x 3 - 3xy 2 ) 


11.41 


11.56 


11.58 


B 2 g 


z{y 3 - 3x 2 y) 


13.34 


13.39 


13.25 


Bin 


y 3 - 3x 2 y 


13.58 


13.58 


13.37 



To evaluate the correlation function k(t), we sample an initial ensemble 
of 1000 walkers from every 100 steps of a VMC walk. This initial ensemble 
is propagated by a DMC sidewalk with a time step of At = 10 Hartree" 1 . 
In the 4 He-benzene system, the energy differences of interest are sufficiently 
large such that we can employ and compare all three DMC implementations 
discussed in Sec. 2.2. For the hybrid branching/weighting scheme, the en- 
semble size and sum of weights in the DMC propagation are kept at approx- 
imately 1000 on average by choosing an appropriate set of DMC parameters 
Wmin,w max , and 77 (Eq. (|l4|)). For DMC with pure weights and DMC with 
pure branching, the only adjustable parameter is the update parameter r\. 
About 500 independent decays k(t) are generated in this manner, and subse- 
quently used as input for the MEM inversion, resulting in the spectral function 
k(E). Each choice of projector results in a single dominant peak in the 
corresponding k(E), and the peak positions are listed in Table |^. These ex- 
cited state energies show general agreement (to within ~ 0.2 cm" 1 ) between 
the three DMC implementations. 

In Fig. |^ we superimpose the spectral functions obtained using the hy- 
brid branching/weighting approach. There, the peaks are grouped in dou- 
blets whose splittings range from ~ 0.2 — 0.7 cm" 1 . These doublets are 
due to projectors which are symmetric and antisymmetric with respect to 
reflection about the benzene plane. They constitute a tunneling splitting 
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E-E Q [cm- 1 ] 

Figure 3. Spectral function k{E) for 4 He-benzene, computed using a hybrid branch- 
ing/weighting approach. Note that this plot represents a superposition of k(E)'s obtained 
from multiple projectors, each yielding a single peak from the MEM inversion. 



between the two equivalent global potential minima along the benzene C%- 
axis, above and below the aromatic ring plane. Tunneling of helium around 
a planar moiety has also been observed in basis set calculations for the 
2,3-dimcthylnaphthalene-He complex, where the magnitude of the splittings 
ranged from < WT 4 cm -1 for localized states up to 3.2 cm" 1 for highly de- 
localized states.E2l The tunneling splittings which we obtain here exhibit a 
decrease in magnitude with increasing energy. Since the energies of highest 
levels correspond to about 12% of the 4 He-benzene potential energy minimum, 
this decrease in the tunneling splitting can be attributed to increasing anhar- 
monicities in the 4 He-benzene interaction potential as these levels approach 
dissociation. Inclusion of the benzene rotational kinetic energy term into the 
Hamiltonian of Eq. (|2l|) qualitatively changes the features of the energy spec- 
trum, removing this decrease in the tunnel splitting. The specific effects of 
this rotational contribution, as well as theugeneral physics of 4 He./v-benzene 
clusters, will be reported in a future study.E!l 

4 Conclusion 

We have extended the applicability of the POITSE method by introducing 
branching processes in the DMC evaluation of an imaginary time correlation 
function k(r). The effects of branching were tested in the determination of 
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excited state energies for two simple systems, namely the one-dimensional 
ammonia inversion mode and the six-dimensional 4 He-benzene van der Waals 
modes. While in an ideal situation one would employ a pure weighting scheme 
in the DMC propagation, in the ammonia study we were faced with the prob- 
lem of evaluating a slowly decaying k(t) corresponding to a small tunneling 
splitting. Thus, the incorporation of branching in the DMC sidewalk is es- 
sential for the stable computation of small energy differences. We have also 
provided a comparison between the various branching schemes and the pure 
weighting scheme in the 4 He-benzene example, and have demonstrated that 
the results obtained are in good numerical agreement. 

The incorporation of branching as described here has been critical in 
allowing excited state energies_ta.now be evaluated for much larger systems 
using the POITSE approach. E3E3 Another current modification in progress 
includes the implementation of descendant weighting techniquesHu'EirEj to 
construct an estimate of the exact ground state wave function |</>o) "on-the- 
fly". The projector A would then operate on the exact \4>o), instead of an 
approximate trial function \^t)- These improvements in the general POITSE 
methodology open the way for efficient and accurate Monte Carlo evaluation 
of excited state energies for large systems. 

Acknowledgments 

Financial and computational support from the National Science Foundation 
through grant CHE-9616615. An allocation of supercomputing time from the 
National Partnership for Advanced Computational Infrastructure (NPACI) is 
gratefully acknowledged. 

References 

1. D. Blume, M. Lewerenz, P. Niyaz, and K. B. Whaley, Phys. Rev. E 55, 
3664 (1997). 

2. D. Blume, M. Lewerenz, and K. B. Whaley, J. Chcm. Phys. 107, 9067 
(1997). 

3. D. Blume, M. Mladenovic, M. Lewerenz, and K. B. Whaley, J. Chem. 
Phys. 110, 5789 (1999). 

4. D. Blume and K. B. Whaley, J. Chem. Phys. 112, 2218 (2000). 

5. R. K. Bryan, Eur. Biophys. J. 18, 165 (1990). 

6. D. Blume, Ph.D. thesis, University of Gottingen, 1998. 

7. B. L. Hammond, W. A. Lester, Jr., and P. J. Reynolds, Monte Carlo 
Methods in Ab Initio Quantum Chemistry (World Scientific, Singapore, 



Copyright © 2002 by World Scientific, reprinted with permission. 15 



1994). 

8. R. Assaraf, M. Caffarel, and A. Khelif, Phys. Rev. E 61, 4566 (2000). 

9. R. N. Barnett, P. J. Reynolds, and W. A. Lester, Jr., J. Comp. Phys. 

96, 258 (1991). 

10. D. Blume, M. Lewerenz, F. Huiskcn, and M. Kaloudis, J. Chcm. Phys. 
105, 8666 (1996). 

11. C. Lcforcsticr, J. Chcm. Phys. 94, 6388 (1991). 

12. A. Nino and C. Munoz-Caro, Computers Chcm. 19, 371 (1995). 

13. C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. Rev. Lett. 60, 
1719 (1988). 

14. M. Snajdr and S. M. Rothstein, J. Chcm. Phys. 112, 4935 (2000). 

15. M. P. Nightingale and V. Melik-Alaverdian, see article in this volume. 

16. Y. Kwon and K. B. Whaley, J. Chcm. Phys. 114, 3163 (2001). 

17. P. Hobza, O. Bludsky, H. L. Selzle, and E. W. Schlag, J. Chcm. Phys. 

97, 335 (1992). 

18. A. Mushinski and M. P. Nightingale, J. Chcm. Phys. 101, 8831 (1994). 

19. K. B. Whaley, Advances in Molecular Vibrations and Collision Dynamics 
(JAI Press Inc., Greenwich, CT, 1998), p. 145. 

20. A. Bach, S. Lcutwyler, D. Sabo, and Z. Bacic, J. Chcm. Phys. 107, 8781 
(1997). 

21. P. Huang and K. B. Whaley, (2001), to be submitted to Phys. Rev. B. 

22. A. Viel and K. B. Whaley, J. Chem. Phys. (2001), submitted to J. Chem. 
Phys. 

23. K. S. Liu, M. H. Kalos, and G. V. Chester, Phys. Rev. A 10, 303 (1974). 

24. J. Casulleras and J. Boronat, Phys. Rev. B 52, 3654 (1995). 

25. M. Hornik, M. Snajdr, and S. M. Rothstein, J. Chem. Phys. 113, 3496 
(2000). 



Copyright © 2002 by World Scientific, reprinted with permission. 16 



